Physiologically based radiopharmacokinetic (PBRPK) modeling to simulate and analyze radiopharmaceutical therapies: studies of non-linearities, multi-bolus injections, and albumin binding

Background We aimed to develop a publicly shared computational physiologically based pharmacokinetic (PBPK) model to reliably simulate and analyze radiopharmaceutical therapies (RPTs), including probing of hot-cold ligand competitions as well as alternative injection scenarios and drug designs, towards optimal therapies. Results To handle the complexity of PBPK models (over 150 differential equations), a scalable modeling notation called the “reaction graph” is introduced, enabling easy inclusion of various interactions. We refer to this as physiologically based radiopharmacokinetic (PBRPK) modeling, fine-tuned specifically for radiopharmaceuticals. As three important applications, we used our PBRPK model to (1) study the effect of competition between hot and cold species on delivered doses to tumors and organs at risk. In addition, (2) we evaluated an alternative paradigm of utilizing multi-bolus injections in RPTs instead of prevalent single injections. Finally, (3) we used PBRPK modeling to study the impact of varying albumin-binding affinities by ligands, and the implications for RPTs. We found that competition between labeled and unlabeled ligands can lead to non-linear relations between injected activity and the delivered dose to a particular organ, in the sense that doubling the injected activity does not necessarily result in a doubled dose delivered to a particular organ (a false intuition from external beam radiotherapy). In addition, we observed that fractionating injections can lead to a higher payload of dose delivery to organs, though not a differential dose delivery to the tumor. By contrast, we found out that increased albumin-binding affinities of the injected ligands can lead to such a differential effect in delivering more doses to tumors, and this can be attributed to several factors that PBRPK modeling allows us to probe. Conclusions Advanced computational PBRPK modeling enables simulation and analysis of a variety of intervention and drug design scenarios, towards more optimal delivery of RPTs. Supplementary Information The online version contains supplementary material available at 10.1186/s41181-023-00236-w.


S.1 Naming Convention for Variables
It is often very common in the community of mathematics to show different variables with single letters such as x, α, etc.However, this approach is not really scalable and becomes very challenging in assigning names to variables when there are lots of them in a model.Hence we developed our own naming convention in developing our PBRPK model which makes it easily scalable.
The naming conventions are quite similar to the strategy when we use objectoriented programming.For instance, suppose that you want to refer to the hot (i.e.labeled) radiopharmaceuticals in the interstitial space of the liver.The name of that variable will be like: Liver.int.H in which "int" reflects the interstitial space, and "H" reflects the hot radiopharmaceuticals.

S.2 Structure of the model
As we discussed in the presented paper, the PBRPK model is created by connecting different organs to each other using ordinary differential equations.Lungs, heart, adipose, skin, bone, red marrow, and muscle are considered as the receptor-negative organs (that do not express any PSMA binding sites), while the liver, spleen, salivary glands, prostate, GI, tumor, and kidney are considered as receptor-positive organs which means that they effectively express PSMA binding sites.
Receptor-negative organs do not have any radiopharmaceutical uptake while the receptor-positive organs have radiopharmaceutical uptakes.This uptake can be in the form of specific uptake (bound to actual receptors), or non-specific uptake (which still can be modeled as binding to some imaginary receptors.The organs in each category have a very similar structure and are only distinguished due to the different parameters that they have.In the following section, we will explain the internal structure of each category of organs.

S.2.1 Internal Structure of the Organs in PBRPK
All receptor-negative organs have a similar structure which is shown in Figure S1.The only exception is the lungs in which the inflow of blood into the organ is from veins (instead of arteries) and the outflow of blood from the lungs is to the arteries.
In contrast, for the receptor-positive organs, we will have slightly different structures to take the receptor binding and internalization process into account (see Figure S2 ) However, due to the complex nature of the kidney, it has its own dedicated structure design to take the kidney functions into account (see Figure S3).Fig. S1: Compartmentalized structure of the receptor-negative organs.Note that there are no "bound" and "intern" compartments, which reflects the fact that the organ is receptor-negative.At each compartment (as the inner structure of the organs), there are generally five distinct species called H, C, AH, AC, Alb which stands for hot and cold (labeled and unlabeled) radiopharmaceuticals, hot-albumin and cold-albumin complex and albumin itself respectively.These species interact with each other which can be summarized in a "reaction graph" (see Figure S4).This reaction graph can be fully described by a system of ODEs as following Fig.S4: Interactions of the different species with each other.The yellow circles represent the reactions between species and those circles labeled with a two-way arrow indicate that the reaction is a two-way reaction.The reactions that are not two-way are the decay process of radiopharmaceuticals.
Using the fact that these species are repeated in each compartment, we can now describe the concentration of each of these species in each compartment by writing a differential equation.Because the ODE for all receptor-negative organs (same as receptor-positive organs) looks similar to each other (except for the numerical values of the parameters as well as the variable names), we explicitly write the ODE system for just one organ from each category here.However, the full details of all of the equations can be found in the table (auto-generated by the Simbiology package in Matlab) attached at the end of this supplementary file (see subsection S.4).
The ODE system for any receptor-negative organ (name it X) will be like this: in which we have Note that the variables inside bracket [•] represent the concentration of that variable.For instance, we have in which we have

S.3 S-Values and Dose Calculation
Understanding the concept of radiation dose and radiation dose rate, as well as their units of measurement, is crucial for calculating radiation dose.Radiation dose, specifically radiation absorbed dose (D), quantifies the total energy absorbed by one unit of mass of a substance.In the SI Unit System, the absorbed dose is measured in gray (Gy), where 1 Gy equals 1 joule (J) of energy absorbed per kilogram of substance.
Commonly used derived units from gray are centigray (cGy = 10 −2 Gy) and milligray (mGy = 10 −3 Gy).The older unit of radiation dose is rad, which represents 100 ergs (1 erg = 10 −7 J) of absorbed energy per gram of tissue or substance.One rad is equivalent to (1/100) Gy or one cGy, and 1 Gy equals 100 rad.On the other hand, radiation dose rate (dD/dt) measures the amount of energy absorbed per unit of time per unit mass of tissue.It can be expressed in units such as mGy (rad) per minute, cGy (rad) per hour, or Gy (rad) per day.
In order to calculate radiation dose, it is necessary to determine the average energy absorbed by 1 gram of tissue in a target organ from the total energy released by the decay of a specific amount of radioactivity.Since x-rays or γ-rays are highly penetrating, a small source emitting these types of radiation, localized anywhere in the body, will irradiate virtually all organs.For instance, if a radionuclide emitting x-rays or γrays is localized solely in the liver (source S), it will deliver radiation dose not only to the liver but also to all other organs (targets T) in the body.Therefore, when a radiopharmaceutical is localized in multiple organs (sources), the radiation dose from each source to each organ (target) must be calculated.Only then can all the contributions be added together to determine the final radiation dose for each organ.
The process of determining radiation dose involves essential steps [1].Firstly, calculate the rate of energy emission (erg/h) from the different types of radiation emitted by the radionuclide within the source volume.Next, determining the rate of energy absorption from these radiations by the target volume.Next, from the above-mentioned, calculating the average dose rate, dD/dt.This dose rate can be captured via the S-Values (shown in Table S1 and Table S2 for Lu-177 self-absorbed doses).Finally, calculating the average dose, D, by multiplying the S-values with the time-integrated activity (TIA).S-values rely on physical data, such as decay characteristics, and organ shape and size.On the other hand, the TIA necessitates radiopharmaceutical distribution data.

Fig. S3 :
Fig. S3: Compartmentalized structure of the kidney.The distinct feature of kidney in filtering the blood leads to its different structure than other receptor-positive organs.

Fig. S5 :
Fig. S5: The structure of PBPK model as implemented in Simbiology Matlab.

Table S1 :
S-Values for tumor and salivary glands.Table adapted from [2].

Table S2 :
[2]alues used for kidney dose calculation.Tableadaptedfrom[2].The following table summarizes the parameter values used for the simulation.Note that the values printed in the "initial values" column are the values used as the initial condition to solve the ODE system.The values printed under the "value" column are intermediate values that are transformed into initial values by the "initial assignments", which are also listed as separate tables in this section.